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We study resonant tunneling through a single-level quantum dot in the presence of strong 
Coulomb repulsion beyond the perturbative regime. The level is either spin-degenerate or can 
be split by a magnetic field. We, furthermore, discuss the influence of a bosonic environment. Using 
a real-time diagrammatic formulation we calculate transition rates, the spectral density and the 
\Q • nonlinear I — V characteristic. The spectral density shows a multiplet of Kondo peaks split by the 

| transport voltage and the boson frequencies, and shifted by the magnetic field. This leads to zero- 

^\ ■ bias anomalies in the differential conductance, which agree well with recent experimental results for 

the electron transport through single-charge traps. Furthermore, we predict that the sign of the 
zero-bias anomaly depends on the level position relative to the Fermi level of the leads. 
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I. INTRODUCTION 



The experimental study of tunneling through zero-dimensional states in quantum dots with high charging eaepsies 
has received recently considerable interestfirLj. Theoretical studies cover the classical regime (high temperatures pV3 as 
well as the quantum-mechanical regime (low temperatures )E3i 2 J. In the latter case, Coulomb blockade and resonant- 
tunneling phenomena together with noncquilibrium generalizations of the Kondo effect are expected to occur. This 
leads to zero-bias anomalies in the differential conductance, which recently have been observed by Ralph & BuhrmanEJ. 
In this article, we present a new real-time diagrammatic approach to describe resonant tunneling at low temperatures 
and compare our results to the latter experiment. 
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■ The experiments for Coulomb blockade phenomena in zero-dimensional systems are usually performed in double- 
barrier resonant-tunneling structuresEm, split-gate quantum-dot devicesotL quantum point-contacts with single-charge 
trap statesE-3 and quite recently also in ultrasmall metallic tunnel junctionsB with Al particles of diameter below lOnm. 
In the latter experiment, the level spacing is of order O.bmeV which is comparable to the Coulomb charging energy in 
usual quantum dots. Therefore, the quantum dot is described by the nonequilibrium Anderson model where the energy 
level e a (with spin label a) is coupled via tunneling barriers to two electron reservoirs with different electrochemical 
potentials fir, and fiR. The charging energy is described by a strong on-site Coulomb repulsion U which suppress* 
double occupancy of the dot level. In equilibrium, it is well known from the theory of strongly correlated fermionsl 
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q that the spectral density of the dot can exhibit a Kondo resonance at the Fermi level. It occurs for a low lying level 
• • i e<7 — fJ-L,R < — r and weak Zeeman splitting \e a — tg\ < Y, where F/^ is th e level width in the noninteracting case, 
. £h \ and provided that temperature is lower than the Kondo temperatureEal 2 ^ Tk = l/2(J/r) 1 / 2 exp [ire(e + U)/(YU)]. 

Since the weight of the equilibrium spectral density at the Fermi level is pjpoortional to the linear conductance, an 
enhancement of the latter due to Kondo-assisted tunneling was predictecOiiS. Typical values for quantum dots are 
U ~ lmeV and Y ~ 50fieV which, for e ~ — Y, yield a Kondo temperature of the order Tk ~ 50mK. Due to heating 
effects such temperatures are still hard to realize in realistic dots. A more pronounced feature was found fo, 
nonlinear conductance which shows a zero-bias maximum even for temperatures above the Kondo temperaturetiaE. 
At zero magnetic field the spectral density of each spin channel exhibits a Kondo resonance at each of the chemical 
potentials. An applied magnetic field causes the Kondo peaks to shift from the chemical potential by the Zeeman 
energy, in opposite directions for the different spin channels. As a consequence, Kondo-assisted tunneling can only 
occux-.if the bias voltage exceeds the Zeeman energy. Therefore, the zero-bias anomaly is split by an applied magnetic 
fieldlij. These features have been observed experimentally by Ralph & BuhrmanEj. They measured the differential 
conductance through single-charge traps in a metallic quantum point-contact. Although this system does not allow 
a controlled variation of the level position, the appearance of a zero-bias maximum with a peak height varying 
logarithmically with temperature clearly demonstrates the mechanism of Kondo-assisted tunneling. However, the 
detailed comparison of the line shape in the experiment and the existing theory showed significant deviations. 

In this paper, we will describe a new approach to calculate nonequilibrium properties of strongly correlated meso- 
scopic systems coupled to fermionic or bosonic baths via particle and energy exchange. It consists in a real-dime 
diagrammatic approach closeljijelated to path-integral methods formulated in connection with dissipationE3~E!j or 
tunneling in metallic junctionstiTLJ. One of the difficulties of the present problem lies in the fact that we have to 
account for the Coulomb interaction in a nonperturbative way. Therefore, the standard many-body diagrammatic 
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approaches are not sufficient (since Wick's theorem cannot be applied naively). We circumvent the problem by keeping 
track explicitly of the time evolution of the density matrix of the dot and tracing out only the bath degrees of freedom, 
which are assumed to be in equilibrium. The final diagrammatic language is set up by an expansion in the coupling 
to the fermionic reservoirs whereas the strong correlations on the local system are exactly taken into account. The 
basic step is the calculation of transition and current rates between different states of the dot. We present an exact 
expression for these rates as the sum over all irreducible diagrams. The transition rate is used to set up a formally 
exact Master equation from which the time-dependent probability distribution for the dot can be calculated. The 
current rate is generally not identical to the transition rate since it contains also the number of particles transferred 
to the reservoir where the current is calculated. This number can take arbitrary values if one considers all higher 
order processes. The occupation probabilities multiplied with the current rates are used to calculate the current 
flowing through the system. In earlier publications |— we have presented this technique in connection with tunnclmc 
through a metallic island with a continuum of statesE-3 and demonstrated the equivalence to path-integral methodsEi 
There we used an approximation for the rate to sum up "inelastic resonant tunneling" processes to arbitrary order 
where different electrons tunnel coherently back and forth between the island and the reservoirs. Here, we apply an 
equivalent approximation to describe resonant tunneling between metallic leads through an ultrasmall quantum dot 
with a single level. An important advantage of our approach is that we can solve the noninteracting limit exactly and 
can control systematically if this limit is contained within a given approximation for the correlated case. The theory 
is current conserving and can be used for the calculation of correlation functions or Green's functions as well. 

For the case where the dot level is M-fold degenerate we recover for M > 2 and a low lying dot level a Kondo 
peak in the spectral function. An applied transport voltage leads to a splitting of the Kondo peak (at /x a where a 
denotes the lead), which results in a zero-bias anomaly in the differential conductance, such that the conductance-has 
a maximum at V — 0. On the other hand, if the dot level lies above the Fermi levels of the leads fi a we predictliij a 
zero-bias anomaly in the conductance which has a minimum at V — 0. 

Several extensions will be considered. We study the case where the (spin-) degeneracy of the dot level is lifted, e.g., 
by an applied magnetic field. In this case the Kondo peaks of both spin channels move apart from each other by the 
level spacing e CT — eg-, and Kondo-assisted tunneling sets in only at transport voltages eV exceeding this splitting. The 
calculated conductance agrees well with the experimental results of Ref. p4[ 

We, furthermore, account for inelastic interactions with bosonic modes coupled to the dot. They describe applied 
time-dependent fields, interaction with phonons, or the fluctuations in the electrodynamic environment. The inves- 
tigation of this field has started only recently in connection with transport through interacting quantum dots. The 
classical regime (r 3> T) has been analyzed for time-dependent fieldscil and bosonic environmentst£l. Photon- and 
boson-assisted tunneling leads here to resonant side peaks in the Coulomb oscillations-svhich can be used to analyze 
the complete excitation spectrum of the dot. The results agree well with experiments^]. In the Kondo regime it has 
been foundEj that time-dependent perturbations split the Kondo resonances which leads to satellite anomalies in the 
differential conductance and offers the possibility to realize pump effects which arc purely based on the presence of 
Kondo resonances. The linear AC-conductance has been analyzed in Ref. [38|. In this paper, we will investigate the 
influence of an external bosonic field on transport phenomena through ultrasmall quantum dots at low temperatures 
and small boson frequencies (compared to T). The emission and absorption of bosons causes additional Kondo sin- 
gularities, for a one-mode field at \x a + nujB, where n = ±1, ±2, . . .. Again, these resonances lead to corresponding 
anomalies in the differential conductance, which are inverted if the level position of the dot is moved through the 
Fermi energy. 



We concentrate here on a dot containing only one energy level e a with degeneracy M . In a magnetic field, due 
to the Zeeman energy the level position is spin dependent. The dot is connected via high tunneling barriers to two 
large noninteracting reservoirs and coupled capacitively to an external gate voltage. The model Hamiltonian of this 
"single-electron transistor" (see Fig. pi) is 



Here, Hp is the Hamiltonian for the dot. It includes the Coulomb interaction of the dot electrons, which is described 
within the capacitance model of a single-electron transistor by the capacitances Cl and Cr for the left and right tunnel 
junction and Cq for the gate. Furthermore, the dot electrons are coupled to bosonic modes u q with electron-boson 
coupling g q . Thus, Hd reads 



II. HAMILTONIAN 
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Y,^ } na + E c (N-n G y + Y,u q dld q + NY,9 q (d q + dl) . (2) 



(Throughout this work, we set % = k = 1 and use the convention e > 0.) The particle number on the dot with 
spin a is n a = c^c a , and N = 'Yl ia n "- The scale of the charging energy is provided by Ec = e 2 /2C, where 
C = Cl + Cr + Cq is the total capacitance of the system. The transistor can be tuned by the gate voltage Vq via 
enc = ClVl + CrVr + GgVg- We remark here, that H is invariant under a global shift of all energies. Therefore, 
we can always choose symmetric bias, Vl — —Vr — V/2. 

The next term in Eq. (|l|), H a = ^2 ka tka^ ac flkaa , describes the reservoir a of noninteracting electrons in the leads. 
Finally, the dot is coupled via tunnel barriers to the left and right lead. This coupling is described by the tunnel 
Hamiltonian 

H T = J2( T >l a c« + h.c). (3) 

kaa 

The bosqniq— .modes can represent interaction with phononsS~0 or fluctuations of the electrodynamic 
environmentc3c3 very similar to the Caldeira-Leggett modeled. For our theory no assumption is needed for the 
specific kind of the modes cu q and the couplings g q . In this way we are able to present a general result for the current 
which shows the influence of inelastic interactions for an arbitrary environment. 

The capacitive model is equivalent to the Anderson Hamiltonian, which we obtain by defining the interaction 
Uq = 2Ec and shifting the level position e^P + 2Ec — cVgCg/C — a c eV/2 — > (V, Vg). The asymmetry factor 
a c = (Cl — Cr)/C accounts for a different capacitive coupling of the left and right lead. We see here, that the effective 
level position in the Anderson model depends on the gate voltage, as well as, due to a c , on the transport voltage. The 
dot is then described by 

H D = J2^\V,V G )n a + U ^ CT '+E^4^+^E^K+4)- ( 4 ) 

c ct<ct' q q 

An unitary transformation0 with V = exp(-iNip) and ip = i T iq (9q/ UJ q)( d \ ~ d q) yields H = VHV- 1 = H Q + H T , 
where H = Hn + Hp , 

Hp = e a n a + U ^ n <y n <y' + ^ w q d l d q ( 5 ) 

a a<a' q 

and 

H T =Y t ^ kaa c a e iv + h.c.). (6) 

kaoi 

The electron-boson interaction renormalizes the level position and the Coulomb repulsion, e CT = ei — ^2 q g q /oj q and 

U = Uo — ^Tliq^ql^qi and the tunneling term acquires phase factors e ±lv . 

For strong Coulomb repulsion U we restrict ourselves to states with N = 0,1. In lowest order perturbation theory 
the rates for tunneling in and out of the dot to reservoir a are 

2^ 7 ±(£) = 2ttJ dE'^(E')P ± (E - E') , (7) 

where 2n^{E) = T a {E)f±{E) is the classical rate without bosons, T a (E) = 2nY /k \ T k\ 2s ( E - e ka), and f+(E) is 
the Fermi distribution of reservoir a with electrochemical potential fi a , while f^(E) = 1 — f£(E). Finally, 

P ± (E) = -!- / dte lEt < e M0) e -M±t) >0 ( 8 ) 



2tt _ 

describes the probability that an electron absorbs (P + ) or emits (P~) the boson energy The exijectation 

value is taken with the free boson Hamiltonian. These probabilities satisfy the condition of detailed balanceea 

P-(E) = P + {-E) = e pE P+{E) (9) 

The classical rates combined with a master equation are sufficient in the perturbative regime r = ^ Q T a <C T.S 
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In order to go beyond the perturbative regime, we need a nonperturbative treatment of the tunneling, where 
quantum fluctuations yield finite life-time broadening and renormalization effects of the dot levels. As an illustration 
we first assume that (for B = 0) the broadening is given by the sum of the classical transition rates, Eq. (0). Using 
Kramers-Kronig we deduce the renormalization and obtain for the self-energy 



a(E) 



dE' 



M 7 +(£ / ) + 7 ~(E v ) 
E-E' + i0+ ' 



(10) 



where 7 = X) Q 7a- The aim of the present letter is to test and extend this simple physical picture within a 
systematic and conserving theory for all Green's functions and the current. To achieve this we use a real-time 
technique developed in Ref. 34 3qJTa which provides a natural generalization of the classical and cotunneling theory 



to the physics of resonant tunneling. 



III. DIAGRAMMATIC TECHNIQUE 



A quantum-statistical expectation value of an operator A at time t is given by 

(A(t)}=tT{p A(t)ff}, 



(11) 



where A{t)fj — exp[iH(t — t )]Aexp[—iH(t — t )] is the operator in Hcisenberg picture with respect to the initial time 
<o- Permutation under the trace yields (A(t)) = tr[p(t)A] with A in Schrodinger picture. The density matrix p(t) 
evolves in time via p(t) = e - 2S ( t -*°V(io)e li?(t ~* o) . We assume that the initial density matrix po — p(to) factorizes 
into parts for the dot electrons, the bosons, and the leads: 



Po 



The leads are treated as large equilibrium reservoirs with fixed electrochemical potentials p c 
describe the electrons in the leads by Fermi functions f a (E) and the density matrix reads 



Pa 



1 



exp [-/3(H a - p a N a )} 



(12) 

-eV a . Therefore, we 
(13) 



where f3 



1/T and N a = ^Zka a kaa a k<ra the number of electrons in the lead a. The normalization factor Z$ is 



determined by trpQ = 1- The boson part reads 



b 1 
Pa = ~zz exp 



Uqd\d q 



(14) 



The temperature of the boson bath, Tb = 1//3b, may differ in real experiments from the electron temperature T. 

For the initial distribution of the dot we assume that it is diagonal in the many-body dot states, \x), which include 
the strong correlations within the quantum dot but are assumed to have fixed occupation numbers, 



PE 



(15) 



with P° = 1. We will see later, that in the stationary limit, i.e., when to is shifted to minus infinity, all the 
physical quantities are independent of the choice of P2. 

In the following, it is convenient to change to the interaction picture with respect to Hq. This implies A(t)g[ = 
Texp (-2 r t tQ dt' H T (<')/) A(t)iTexp (-i Jl dt' H T {*)i) 

in which T is the time-ordering and T denotes the anti- 
time-ordering operator. We write the integrals as one contour integral J„ dt 1 . . . over the Keldysh contour. It is 
parameterized by the "time" t' which first runs forward from to to t and then backward from t to to- In the 
diagrammatic language, the Keldysh contour is represented by horizontal lines running from the left to the right and 
then back to the left (see Fig. ||). We find 



(A(t)) 



tr 



po T K exp 



dt'Hrit')! A(t) 



K 



(16) 
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Here, we have introduced the Keldysh time-ordering operator Tr, which orders all following operators along the 
Keldysh contour such that the one with the later "time" along the Keldysh contour appears at a more left position 
(without any sign change for an exchange of Fermi operators). 

In the following we will encounter also higher order correlation functions of the type (Tk A\(tx)iA2{t2)i ■ ■ ■ A n (t n )j). 
The final time t of the Keldysh contour is then given by max{ti, . . . , t n }. 

For a diagrammatic description we expand the exponential with respect to the tunneling Hamiltonian and obtain 

f>i) m / dt[[ dt' 2 ...[ dt' m T K {H T (t[) I H T (t' 2 ) I ...H T (t'J I f[A l (t l ) I \ 

Jk Jk Jk { ti J J (17) 

t[ > t' 2 > ... >t' m 

in which the relation t' l > t' 2 > . . . > t' m has to be understood with respect to the Keldysh contour. The time-ordering 
operator Tk acts also on the operators Ai(ti)i and puts them on the right place between the tunneling Hamiltonians. 

The next task is to perform the trace of each term of the expansion. We insert the tunnel Hamiltonian (^) and notice 
that the Hamiltonian Hq is bilinear in the lead electron operators. For this reason, Wick's theorem holds for these 
degrees of freedom, i.e., the lead electron operators are contracted in pairs. This includes contractions between pairs 
of field operators from Hj> as well as contractions to lead electron operators which can be present in the operators Ai . 
These contractions are given by equilibrium distribution functions. For the dot electrons, the situation is different. 
The Coulomb interaction is expressed by a quartic term of dot electron operators. Therefore, Wick's theorem does 
not hold for this part of the system. A product of dot electron operators can not be contracted into pairs, but has 
to be treated explicitly. The trace over the bosonic part, however, poses no problem. Each tunneling term contains 
an exponential e ±JV of the bosonic operators (p. The trace over the product of such exponentials can be performed 
easily (see Eq. (p^)) provided that the operator A also contains only such exponentials. 



71 

(T K Y[MU)) =tr 

z=l 



A. The rules 



With regard to the applications discussed in section 1IIB and 1IIC, we assume here that the operators A, 
senting "external vertices" ) depend on the lead and boson degrees of freedom only in the form 



Ai — Ai 




(repre- 



(18) 



Each term of the expansion Eq. ( |l7| ) is visualized by a diagram (see Fig. g). There is a forward and a backward 
propagator symbolized by the upper and lower horizontal line, running from i to t and back from t to t , respectively. 
Along this time path, we arrange internal and external vertices according to their time-ordering. The internal vertices 
emerge from the insertion of the tunneling Hamiltonian (^) into the expansion ((l?]). Each of them corresponds to 
a product of a lead and a dot electron operator and a phase factor e ±lv . After integrating out the lead degrees of 
freedom, all vertices (either internal or external) containing a lead electron operator are connected in pairs by directed 
tunneling lines (dashed lines) j%(t, t') from t' to t, with *f*(t, t') = j+(t - 1') for t < t' and jg(t, t') = 7 ~(t - t 1 ) for 
t > t' with respect to the Keldysh contour with 7„(t) = / dEe~ lEt ^{E). These tunneling lines represent contractions 
of lead electron operators. 

There are vertices from which a tunneling line leaves (representing a\ lja (t)c a (t)e lip ^ which removes a dot electron 
with spin a) and others to which a tunneling line enters (visualizing c] J {i)ak a a{t)e~ lLp ^ which adds a dot electron 
with spin a). Fermi statistics, furthermore, yield a minus sign for each crossing of tunneling lines. 

In the interaction picture the dot electron operators get exponential factors which contain the energies e x of the 
many-body dot states x given by e x \x) = Hd\x)- The order of the electron operators may induce furthermore a 
minus sign due to Fermi statistics. 

The trace over the boson operators gives rise to a factor of the form 



Cb (tl,t2, ■ ■ ■ , t m , ti , t 2 ) • 



.,t' m ,) = {T K e-MVe 



itp(t 2 ) 



-i<p(tm) p i<p(t'i) p i<p(t' 2 ) 



Since if is linear in the boson operators, we get 



(19) 



(20) 



l<3 



1<J 



1,3 
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We write P K {t,t') = P + (t,t') for t < t' and P K {t,t') = P~(t,i') for t > t' on the Keldysh contour with P ± (t) = 
J dEe~ lEt P ± (E). In the diagrammatic language, we represent the factors P K by boson lines connecting each vertex 
with each other. 

A summary of these rules are given in Appendix |a|. 

In order to calculate stationary transport properties it is convenient to change to an energy representation. Without 
loss of generality we assume that the times ti, . . . ,t n of the correlation function ( |l7|) are ordered on the real axis 
according to t n < t n -\ < . . . < ti = t. This may be different to the ordering on the Keldysh contour which depends 
on whether the times lie on the upper or lower branch. In the stationary limit we can set to — — oo and t = ti = 0. 

We consider the Laplace transform 

G(E 2 ,E 3} ...,E n ) = H)"- 1 f dt 2 f 2 dt 3 ... f"" dt n e lE ^e lE ^ . . . e lE ^ (T K A 1 (0)A 2 (t 2 ) . . . A n (t n )) . (21) 

J— oo J— oo J— oo 

We will account for the exponential factors exp(i£"iti) (i = 2, . . . , n) by drawing directed virtual lines from the external 
vertices with time ti to the last vertex with time t\ = and assigning the energy Ei to this virtual line. 

Performing the time integrals we end up with diagrammatic rules in energy representation. These rules are sum- 
marized in Appendix ||. 



B. Master equation and stationary probabilities 

In this section we will derive a formally exact expression for the central object of this paper: the quantum-mechanical 
transition rate £ X ' iX (£', t) for the reduced system (the dot) to go from a state x' at ti me t' to a state \ at time t. This 
rate will serve as an input for a formally exact and time-dependent Master equation which, in principle, could be used 
to calculate all occupation probabilities of the dot as a function of time for an arbitrary initial state. Similar Master 
equations are welUuiown and successfully applied in connection with macroscopic quantum coherence phenomena in 
spin bose modelsETEl The connection to these path-integral approaches will be described in Ref. |5f| 

A matrix element of the reduced density matrix of the dot at time t, P^{t), is given by the quantum-statistical 
expectation value of the projector (|x2)(xi|)M 

P£(t) = {(\X2)(xi\)(t)), (22) 

i.e., we have to set n = 1 and A\ = |x2)(xi| in Eq- (fl7|). The reduced density matrix commutes with the particle 
number on the dot. Thus, the operator |x2)(xi| is unaffected by the unitary transformation and no tunneling and 
boson line emerges from this external vertex. The matrix element P* 1 (t) can be expressed by the reduced propagator 

II*, 1 '** (t', i) from x'i at time t' forward to xi at time t and then from X2 at time t backward to x' 2 at time t' 

E^(On&£ (23) 

Xi 5X2 

The propagator is the sum of all diagrams with the given states at the ends and can be expressed by an irreducible 
self-energy part ^y^* {f, t), defined as the sum of all diagrams in which any vertical line cutting through them crosses 
at least one tunneling or boson line. The propagators for the four lines attached to the self-energy are not included 
in T, x }' Xl (£', t). We obtain an iteration in the style of a Dyson equation (see Fig. 51), 

X2'X2 I — 1 

ng;S(*',i) = n(°)^( i ' > t)^, x ^ X2 ^ + J* dt * fj d *i u Hi (*'' *0 E S';S fa' n(0) £ (*2» *) . (24) 

where n(°>£(f,i) = exp[— i(e Xl - e X2 )(t - 1')] is the propagator of the isolated quantum dot. Multiplying this 

equation with P x Ht'), summing over the states x'i;X2 and differentiating with respect to t, we obtain together with 

X2 

Eq. (E3T) and setting t' = to 

_ .1 _ 1 ** tn 



d 
~dt 
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This formally exact equation is the most general kinetic equation for the reduced density matrix of the dot. No 
assumption is necessary for the initial state and the integral kernel £ on the r.h.s. shows that memory effects are 
fully taken into account. 

The equation simplifies considerably if we assume that the initial density matrix is diagonal. In the general case, 
this does not imply that the reduced density matrix stays diagonal for all times. However, for the special case of the 
Anderson model considered here, spin conservation implies that the reduced density matrix will be diagonal for all 
times t > to. Hence, we consider E X ', X = ^-' x anc ^ °btain from (|2^ 



d 
dt 



../ Jtn 



where P x (t) = -P? W denotes the probability to be in the state \ a t time t. For a time-translational invariant system, 
the time-dependent rates £ x ' lX (t', t) depend only on the time difference £ x ' lX (t' — t). Performing the Laplace transform 
of Eq. ( p6| ) one can then study the time evolution of arbitrary initial probability distributions into the stationary state. 

By attaching the rightmost vertex of each diagram E to the upper and lower propagator the minus sign for each 
vertex on the backward propagator yields E x ' iX (i', t) = 0, which allows us to rewrite Eq. ( p6| ) in the form 

jP x {t) = E f t dt ' iPx>(t')Vx>,x(t',t)-Px(t')Z x , x ,(t',t)} , (27) 



x'#x 



We obtain the structure of a master equation with transition rates given by £ x / jX (i', t). 
The stationary distribution is given by 



= Jim P x (t)= lim P x (0) (28) 



tn — > — oo 



and is not the equilibrium one if the electrochemical potentials of the leads arc different. From (26) and ( |27| ) wc 
obtain 

= E P x'^x',x = E t P x' S *',x - P x^x.,x'} , (29) 
x' x'+x 

where 

£ x ,, x = z f dt'X x ,, x (t',0) (30) 

J —oo 

can be calculated directly by using our diagrammatic rules in energy space. Th e prefactor i in (BC|) together with the 



771—1 time integrations over the internal vertices in E gives a factor i m in (Bl) which cancels the factor (— i) m from 
rule 5. 

The self-energy part E x / jX is purely imaginary. This can be seen by changing the vertical positions of all vertices 
on the Keldysh contour (without changing their horizontal position) and reversing the direction of all tunneling and 
boson lines. Consequently, only the energy differences AEj from rule 2' will change sign. Since the number of vertices 
is even (all tunneling vertices are coupled in pairs by tunneling lines), there is no sign change due to rule 5' and the 
number of resolvents is odd. Thus, the whole diagram has been changed to its conjugate complex up to a sign. 

C. The tunneling current 

The tunneling current flowing into reservoir a is defined by I a (t) = e-^(N a (t)) — ie{[H, N a ](t)), which is equivalent 

to 

I a (t) = -*eE{ T * <(«La^O(t)> - rr<(4a fc<ra e-^)(i)>} ■ (31) 

ka 

The tunneling current is an expectation value of a product of a dot, boson and reservoir electron operator (see 
Fig. |). We obtain 
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/«(*) = *£ f dt'p x ,{t')^ x {t',t) = -eY J f dt'p x ,(t')^ x {t',t) 



(32) 



where the partial self-energies E"^ (t', £) are parts of the total self-energy 

z X 'At',t) = Y,{^J t '>^ + K'J t '^} • ( 33 ) 

a 

They describe processes in which the rightmost tunneling line corresponds to reservoir a and is an outgoing (incoming) 
line if the rightmost vertex lies on the upper propagator or an incoming (outgoing) line if the rightmost vertex lies 
on the lower propagator. Their physical meaning is displayed by the current formula ( |32"|) which shows that they give 
the total contribution to the current rate. We can relate them to an intuitively more physical object, namely the 
rate (i', i), p = 0, ±1, ±2, . . ., which describes the transition rate where p particles are transferred to reservoir a. 
Within our graphical language S^f x (t', t) is given by all diagrams where the number of tunneling lines with reservoir 
index a running from the forward to the backward propagator minus the number of tunneling lines with reservoir 
index a running from the backward to the forward propagator is given by p. We obtain 

E^^-^iEE^^c' 4 )' ( 34 ) 

X XV 

This relation together with current conservation is proven in Appendix The factor p shows clearly that E Q± 
describes the contribution to the current rate. In contrast to lowest order processes, i.e., the golden rule rate, where 
p can only take the values ±1, p can be arbitrary for higher order processes. Nevertheless, Eq. ( |33| ) shows that the 
current rate can be calculated as a partial selection of diagrams already contained in the total transition rate £ x ', x - 

We emphasize that the current formula (]3^) together with the Master equation (|2^) constitutes a complete theory 
to describe time-dependent phenomena starting from an arbitrary diagonal initial state. The original problem has now 
been shifted to the evaluation of the various self-energy diagrams which correspond to transition and current rates. 
The self-energies are defined by a set of irreducible diagrams and thus their corresponding perturbation expansion in 
the number of tunneling lines is a well-defined series and contains no divergent time-integrals. 

For time-translational invariant systems the current rates E"^ (t',t) depend only on the time difference t' — t. To 



calculate the stationary current we define in analogy to (p0| 



df^Jt',0) (35) 



x',x / x',x x 

J —oo 

which again can be calculated directly with our diagrammatic rules in energy space. The stationary current is then 
given by 

# = — E P xK'x = ie H P ^x~x ■ ( 36 ) 



D. Green's functions 

After the unitary transformation the Green's functions of the dot electrons read 

G>(t,t') = -i((c„e i *)(t)(cU- i ' p )(t')) (37) 
G<(t,t')=i((cU- iv )(t')(c*e i '?)(t)). (38) 

Here G> and G< are independent quantities since we do not assume equilibrium. For time-translational invariant 
systems, the Green's functions depend only on the time difference G(t,t') = G(t — t'). The Fourier transform 
G(E) = j dte lEt G{t) can be written in the form 

G>(E)=2ilm(-i) [ dte- iEt (T K (c„e i v)(0)(cle- i '?)(t+)) (39) 

J —OO 

G<(E) = -2tlm(-t) f dte- lEt {T K {c a e tv ){Q){cle^)(t-)) , (40) 

J —OO 
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where t means that the the time t lies on the upper (lower) branch of the Keldysh contour. Note that the time 
ordering is defined here by a pure ordering along the Keldysh contour without any sign change if we interchange 
fermion operators. The integrals can be calculated like Eq. ( pl| ) whereby instead of assigning the energy — E to the 
virtual line connecting the external vertices one can change the direction of the line and assign the energy E. 

In order to relate the current to the Green's functions of the dot we consider the first diagram on the r.h.s. of 
Fig. |] (the second one is just the conjugate complex). The external vertex can be either contracted by a tunneling 
line to the upper or lower propagator and we recover immediately the structure. ai the Green's functions G > and G < , 
respectively (see Fig. H). We recover for the stationary current the relationEjiijo 

I s J = -ieJ2 [ dE{ 1 +(E)G>(E)+ 1 -(E)G<(E)}. (41) 

In the case that the couplings to the leads have the same energy dependence, T a (E)/T a i (E) — A QiQ ,<, this can be 
written in the form (which was already derived in Ref. |l4| ), 

7 " - e EE f dE ^^'^ pAEmiE) - /+(£)] ■ (42) 

Q< _Q> 

Here, we used the relation between the Green's functions G^ , G> and spectral density p a = " 27ti " . 



IV. RESULTS 



What we have done so far is to derive a diagrammatic language which allows a systematic description of transport 
processes. We have, furthermore, shown how the physical quantities of interest, the stationary probability distribution 
and the current, can be obtained if we know the value of special diagrams. In this section, now, we will explicitly 
calculate the value of the corresponding diagrams. 

We consider here the case of strong Coulomb repulsion U, i.e., we restrict ourselves to the states with N = 0, 1. 
Diagrams in which a higher occupancy occurs do not contribute since they have resolvents of the order 1/U. 

In the following, the index a labels the singly occupied state with spin a — 1, . . . M. The label \ additionally allows 
an empty dot, x = 0, 1, . . . M. 

In general, we can not sum up all possible diagrams. Therefore, we have to find a systematic criterion which 
diagrams should be retained and summed. 

The simplest approximation is to neglect all diagrams where two or more tunneling lines overlap in time (sec the 
leftmost diagram parts in Fig. |J). This means we include those processes which are also described by the master 
equation with rates obtained in lowest order perturbation theory (sequential tunneling), which is a good description 
at high temperature, T <^T. 

In situations when sequential tunneling is suppressed by Coulomb blockade, the lowest order contribution to the 
current arises due to cotunncling. The rates for a process in which an electron enters the dot from the left lead and 
leaves to right one is described by diagrams with two overlapping lines (see the diagram part in the middle of Fig. |^). 

At lower temperature the perturbative approach is not sufficient. Higher-order processes become important. In 
generalization to cotunncling we have to take into account irreducible diagrams with an arbitrary number of correlated 
tunneling processes, i.e., we include resoriattLtunneling. 

Similar to the case of metallic islandsaEj we proceed in a conserving approximation, taking into account non- 
diagonal matrix elements of the total density matrix up to the difference of one electron-hole pair excitation in the 
leads. The graphical representation of this constriction is that only diagrams in which any vertical line will cut at 
most two tunneling lines are taken into account. 

We give two arguments why this class of diagrams is the most important one. Firstly, since we treat the leads 
as large equilibrium reservoirs there should be a tendency of the system to stay close to diagonal states. Secondly, 
our approximation contains the exact solution for the noninteracting limit U = 0: if there is no electron-electron 
interaction in the dot, electrons with different spin do not influence each other, so that this limit is described within 
our model by choosing M = 1. In this case, the selected diagrams are the only contributing ones. The sum of all 
other, more complicated, diagrams is zero. 

Furthermore, we include only boson lines between vertices which are already connected by tunneling lines, i.e., 

m 

C B (t u t 2 ,...,t m ,t' 1 ,t' 2 ,...,t'J^l[P K {t i ,t' i }. (43) 

2 = 1 
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where the pairs i,-,t^ are already coupled by tunneling lines running from t[ to ti. This amounts to a dressing of the 
tunneling lines 7—5-7. This approximation, while neglecting many diagrams, describes well the spectral density of the 
dot at resonance points. The reason is that position and value of the peaks of the spectral density are determined by 
a self-energy a (see Eq. (ff9l)) which is calculated here in lowest order perturbation theory in V including the bosons. 
Higher orders are small for high tunnel barriers. 

First, we relate the rate S x '. x to an irreducible diagram labeled by </>*r*' (a, <r, (see Fig. |). It has an open 
tunneling line associated with tunneling of an electron with spin a in the junction a carrying the energy E. The 
line is directed from the right to the left and its value together with the corresponding resolvent is included in <j). 
The self-energy is then constructed by attaching the open tunneling line of these diagrams to the upper and lower 
propagator (see Fig. [?]) with the result 



2zlm f dB^2{( X |^|xi>^H«.^^)-0ci|c«r|x>^ 1 («.^^)} 

C,Q XI 

= E{ S ^x + ^x}- ( 44 ) 

a 

where the current rates S"^ correspond to the first and second term, respectively. Again we have made use of the 
fact that a diagram becomes the conjugate complex if we change the vertical position of all vertices and the direction 
of all tunneling and boson lines. As pointed out in Appendix any approximation for cf> will lead to a current 
conserving theory. 

We construct the diagram <j> by iteration (see Figs. || and ||). To do so, we need the diagram ir(E) which is the 
propagator ir(E) while a tunneling line with energy E is running in parallel from the right to the left. This diagram 
can also be expressed as an iteration in the style of a Dyson equation (see Fig. [l(]) 

Xi >x 2 

In analogy to E, the self-energy a(E) denotes the sum of all irreducible diagrams with a tunneling line going 
backward in time. Here, the free propagator in energy space is given by 



^(°) X1 (E) = (46) 

E-(e Xl -e X2 )+iO+ 



Hence, we can solve Eq. ([45|) and find in matrix notation the general relation 

n(E) = [[n^(E)}- 1 -a(E)}- 1 . (47) 

Because of the restriction to two charge states, only the matrix elements Tr a (E) = ttq'q(E) of ir(E) and cT a (E) = 
CTq'q (E) of cr(E) are involved, and we deduce from Eq. ( [47| ) 

n a (E) = — . (48) 

Since at most two tunneling lines are allowed at once, the irreducible self-energy a a (E) consists of only one tunneling 
line. We calculate all contributions, which are depicted in Fig. [Ll[ and get 

J E — E' + i0 + J E - E> + e a > - e a + i0+ v ' 



In the spin degenerate case, this is exactly the relation ( |10[ ) which we found from intuitive arguments. 
According to our rules, Figs. and lead to the self-consistent equation for the diagram (f>(a, a, E) 



4>%(a,a,E)=n°(E) 



it 



^(E)-7a(E)Yl J dE '^-WT^^ {a '^ El) 

a.' 

^EE / dE ' E-E' + e!,-e a+ iO + ^( a '°'^ 



(50) 
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and 



E-E 1 



■iO- 



<P a', a (« ^ > E ) 



(51) 



The stationary probabilities and the current are derived from Eqs. (|29|) and (|36|). To calculate the rates we specify 
Eq. (U) and obtain 



E«+ = 2ilm^ f dE<t>^(a,a,E) , 



(52) 
( 53 ) 

whereas all other rates are zero. 

The correlation functions can be calculated from the diagrams shown in Figs. |l2] and |l3|. We have to consider only 
the latest (i.e. rightmost) correlated part of the diagram. The processes before end up with probability P^ in a 
diagonal state \- We have used the same criterion as for the calculation of the density matrix with one exception. If 
a vertical line lies between the external vertices we allow a cut through at most one tunneling line. Here we have used 
the fact that such a vertical line will in addition always cut the virtual line connecting the external vertices. The sum 
of all these diagrams gives (where we can combine always two diagrams to the imaginary part of one of them) 



G>(E) = / dE 1 G>{E')P+{E' -E) 



G<{E)= dE'G<{E')p-{E' -E) 



with 



G>(E) = 2ilm{ir a (E) 



p o st -EE/ dE 



, P «V*°£ (a, a', E>) + £ ct „ P^,K4 («' CT '< E ' 



E-E 1 



i0+ 



G<(E) = -2ilm{ ir a (E) 



P? + Y, / dE 



, P st <0 («, E ') + Ea> PpK'flfa ^ E ') 



E-E' + i0^ 



(54) 
(55) 

(56) 
(57) 



In the following, we discuss for transparency the effect of the coupling to bosons and the presence of a magnetic 
field separately. 



A. Boson-assisted tunneling 



For zero magnetic field, i.e. e a = e for all ex, we can perform the resummation of the corresponding diagrams for 
the rates and the Green's functions analytically (details are given in Appendix |5]) and find 

FJ = 2neMj2 J dE fr-(E) 7 +(E) - 7^)7+ (£)] \t:(E)\ 2 (58) 

a' 

with tt(E) — 7r a (E). We can write this equation in a more intuitive way by inserting the definition (^) for 7^ 

7 « = lY,J dE J dE' {T a ,, a (E\E)f a ,(E')(l - f a (E)) - T a , a ,(E,E')f a (E)(l - f a ,(E'))} , (59) 

a/ 

where 

T a , a ,{E,E , )=MT a {E)T ol ,{E') f dE 1 P+(E 1 - E)P~(E 1 - P>(Pi)| 2 (60) 
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can be interpreted as a transmission probability for an electron to start from reservoir a with energy E and end in 
reservoir a' with energy E' . From the detailed balance condition (^) we get 

T a , ta (E', E) = e^ E '- E ^ a ,(E, E') . (61) 

This guarantees that the current is zero if all chemical potentials of the reservoirs are identical. 

However, the interpretation of T a ^ a i as a one-particle transmission probability in analogy to a generalization of the 
Landauer-Buttiker formula to inelastic interactionsEil is not correct. We see that the transmission probability still 
depends on the Fermi distribution functions via the self-energy a(E) in the denominator of the propagator tt(E). This 
reflects the many-particle aspect of the electron-electron and electron-bospa interaction in our model. 

Comparing our result for T QiCe / with other approaches in the case M = iEJeJe3, we see that the energy dependence 
of <j(E) has been neglected in all previous treatments. We find that even in the M = 1 case, the energy dependence 
of cr(E) cannot been neglected if the temperature T and the typical frequency ojb of the bosons are smaller than T. 

Without bosons, the current formula is exact up to order 0(T 2 ), i.e., sequential and electron cotunneling are fully 
taken into account. With bosons, cotunneling is not described correctly since we have treated the bosons only by a 
dressing of the tunneling lines. This means that our approximation is not valid in regions where the current is very 
small. However, at resonance we believe our treatment to be correct since there we expect that sequential tunneling 
will be just modified by a renormalization and broadening of the local state of the dot which is described by the 
self-energy u{E) which is calculated in lowest order in V here. Higher orders will be small for high tunneling barriers. 

Finally, we calculate the Green's functions and find 

G > {E) = -2m J dE'-f-(E')P+(E' - E)\tt(E')\ 2 (62) 

G< (E) = 2m J dE' 7 + {E')P- (E' - E) \v{E')\ 2 . (63) 

In equilibrium, i.e. fi a — for all a, we obtain the correct sum rule G > (E) = — exp (j3E)G < (E). Furthermore, for 
the M = 1 case, particle-hole symmetry is satisfied. The spectral density has the form 

p(E) = J dE' [j+(E')p-(E' -E)+ 1 -{E')P+{E' -E)] \tt{E')\ 2 . (64) 

The effect of the resonant-tunneling processes is described by the resolvent ir(E) containing the self-energy <r(E), 
Eq. (p9|). The real and imaginary part of the self-energy express energy renormalization and broadening and determine, 
therefore, the position and the width of the maxima in the spectral density. 

To proceed we consider from now on a one- mode environnaent-|(Einstein model) with boson frequency oj q = lub- 
Experimeptally realizations of this model are optical phononsE311J or by fluctuations of an external LC-circuit with 
frequencyc3~E3 lob = (LC)~ X I 2 . The results for a general environment can be anticipated approximately from the 
one-mode case by a superposition. Furthermore, we choose the special case of two reservoirs a = L/R and constant 
level broadening T/2 = Tl = Tr. 

Defining g = J2 q 9q/u B we obtain P ± (E) = J2nPn$( E =•= ulu b), where 
p n = e~ 9 ^ 1+2Na ^ B ^e nuJB/2TB I n (2gN (ujB)e UJB/2TB ) is the probability for the emission of n bosons with frequency 
ujb- Here, Nq(ujb) is the Bose function and /„ the modified Bessel function. Using Eq. (|49|) we obtaintj 



n..» 1 r,.j:(^-^[>»(^)-i«..(i + *^±^'| 



lma(E) = -Tr^2 Pn [Mj + (E + nLU B )+T(E-nujB)} ■ (66) 

n 

Here ^ denotes the digamma function, and we have chosen in the energy integrals a Lorentzian cut-off at Ec- 

The real part of cr(E) renormalizes the level position to higher energies. Furthermore, it depends logarithmically on 
energy, temperature, voltage and frequency. These logarithmic terms are typical for the occurrence of Kondo peaks. 
Hence, we anticipate logarithmic singularities either for M > 2 or for p n ^ p~ n . This includes not only the degenerate 
case but also the case of a single dot level without spin (M = 1) since the probabilities for absorption and emission of 
bosons are different. It is important to remark here, that for systems coupled to classical time-dependent fieldsEj the 
situation is different since then both probabilities are equal. At low enough temperatures we obtain logarithmic peaks 
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in <t(E) at E = fi a + tiujb (n / for M = 1). They lead to maxima of the resolvent ir(E) at E — /i a + tilob {n > 
for M = 1, n > for M > 1) for e < and at E = fi a + nu B {n < 0) for e > 0. The spectral density d&j] ) shows 
resonances at the same points but, due to the additional functions in the integrand, they are shifted by multiples 
of u>b- This boson- assisted tunneling is completely independent from the influence of the bosons on the self-energy 
a{E). 

The spectral density at different voltages for a low lying level e < is depicted in Fig. [14]. Without an applied 
bias voltage, we obtain (for M = 2) the usual Kondo peak near the Fermi level (which we choose as zero energy). 
The emission of bosons leads to additional resonances at multiples of ll>b- For M = 1 and e < 0, resonances occur 
for negative energies and in the case e > we find resonances at positive energies. In these cases, the effects are 
less pronounced and are only visible for very low temperatures. At finite bias voltages all peaks split and decrease in 
magnitude. 

The resonances in the spectral density can be probed by the nonlinear differential conductance as function of the 
bias voltage V, as shown in Fig. |l5| for the case e < 0. The splitting of the Kondo peak leads to an overall decrease 
of the spectral .density in the energy range \E\ < eV (see inset of Fig. [l5|). For this reason, the conductance shows 
the well-knownoE3o maximum at zero bias. The emission of bosons produces a set of symmetric satellite maxima. 
They can be traced back to the fact that pairs of Kondo peaks can merge if the bias voltage is a multiple of the boson 
frequency (see Fig. |lj). This gives rise to pronounced Kondo peaks at E — ±eV/2 and thus to an increase of the 
spectral density with bias voltage near these points. 

Fig. |l6| shows the differential conductance for e > with and without bosons. A striking result is that the 
whole structure is inverted compared to the case e < 0, and we find a zero-bias anomaly although the Kondo 
peak at zero energy is absent. The coupling to bosons yields satellite steps at \eV\ = ulob- The contributions of 
sequential and cotunneling lead, compared to resonant tunneling, only to a weak bias voltage dependence of the 
differential conductance. This shows clearly that the influence of the logarithmic terms in a{E) are still important. 
The logarithmic peaks in Re u{E) decrease with increasing bias voltage and approach the value of E — e if e is large 
enough. Thus the value of E — e — Kea(E) decreases which in turn leads to an overall increase of the resolvent ir(E) 
and the spectral density p(E) near zero energy (see the inset of Fig. |l^). ■— ■ 

Zero-bias minima- are known from Kondo scattering from magnetic impurities £3 They have been observed, ia. 
recent experiments^ and have been interpreted as 2-channcl Kondo scattering from atomic tunneling systcmsc3'E£l 
or by tunneling into a disordered metaLED Here we have shown that zero-bias minima can also arise due to resonant 
tunneling via local impurities if the level position is high enough such that we are in the mixed valence regime. 

Finally, we have investigated the differential conductance at fixed bias voltage as function of the position of the dot 
level, which experimentally can be varied by a gate voltage coupled capacitively to the dot (see Fig. [l7|). The result 
shows a (classical) pair of peaks at |e| = eV/2 together with satellites (due to emission and absorption of bosons) and 
peaks for |e| > eV/2 (only due to absorption). The energy dependence of Ima(E) gives rise to an asymmetry of the 
peak heights. The peak at e = eV/2 is higher than the one at e = —eV/2 since |Imcr(i?)| = tt\A1^ + (E) + j~(E)\ is 
smaller for higher energies (except for M = 1 when particle- hole symmetry holds). This significant effect is due to 
the broadening of the spectral density by quantum fluctuations. 



B. Magnetic field dependence 



In this section we discuss the effect of an applied magnetic field and do not take into account the coupling to bosons. 
Again we consider the case of two reservoirs and constant level broadenings. Since the energy levels e CT are now spin 
dependent, we can no longer solve the self-consistent equations analytically but have to solve them numerically. 

We find Kondo resonances in the spectral density p a (E) at energies E = \i a + e a > — e CT with a' ^ a. This is due 
to the fact, that the correlation functions G<(E) and G>(E) are mainly determined by the resolvent TT a (E) (see 
Eqs. ( |56| ) and ((5^)) which contains via the self-energy a a (E) logarithmic singularities at the corresponding energies, 



( E > 



c 



Rea°(E) = Y-^Y [in, 

y ' ^ 2n ^ L V 2ttT 



, 1 E 
- Rt<y< | - +i — 



2^T 



Ima a (E) = -7T 
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(E)+Y,J + (E- 



(67) 



(68) 



From Eq. (|42j) we see that only energies within the window defined by the difference of the Fermi functions contribute 
to the current. For this reason, there is no Kondo-assisted tunneling at low transport voltage but sets on if transport 
voltage and level splitting are equal. Therefore, for low lying levels the conductance weak at zero bias found in the 
previous section now splits up into two peaks separated by the twice the level splitting)!-!! (see Fig. [l8|). 
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Ralph-and Buhrman recently measured Kondo-assisted tunneling via a single charge trap of a point contact tunnel 
barrier Ej We follow the model proposed by the authors interpretating the experiment as a realization of the Anderson 
model with strong Coulomb repulsion such that double occupancy does not occur. However we think that the 
interaction energy U and not the conduction band width is the relevant cut-off in this situation. 



A comparison of the experiment and our theory is given in Fig. 18, Fig. [19| and Fig. |20|. We find a good agreement 
for the peaks induced by Kondo-assisted tunneling processes if we set the cut-off U — 30meV. The authors suspect 
the single-charge trap to be a dangling bond, for which they expect U = lOOmeV. Our result agrees in the order of 
magnitude, it gives a hint however, that the state may have a larger extension than an ordinary dangling bond or 
that there is screening due to the copper electrodes or both. The peaks for larger magnetic fields show, however, a 
stronger broadening than predicted from our calculation. Nevertheless, our theory reproduces the experimental curves 
much better than the fits given in Ref. H using perturbation theory since we have taken into account nonperturbative 
effects which are obviously important here. 

The model proposed by the authors of Ref. [24] explains the broad peaks at large voltages by the matching of the 
energies of the empty and the singly occupied dot. Our calculations for this case, however, lead to a broader and 
lower peak for positive voltages in comparison with experiment (see Fig. |l9|). We think, therefore, that due to the 
capacitance asymmetry the system becomes doubly occupied before the empty state is energetically favorable. The 
capacitance asymmetry a c makes then the corresponding resonance peak sharper. An energy dependent transparency 
of the barriers could then explain the different heights. A generalization of our theory to situations, where multiple 
occupancy of the dot is important, is currently under way and will be presented elsewhere. 

Finally, we consider the case when the energy level is above the Fermi energies of the leads. The zero-bias minimum 
found in the previous section splits for finite magnetic field into two minima separated by twice the level splitting (see 
Fig.pl]). 



V. CONCLUSIONS 



In conclusion, we have studied low-temperature transport in the nonequilibrium Anderson model with bosonic 
interactions. The latter yield new Kondo resonances in the spectral density which can be probed by the measurement 
of the nonlinear differential conductance. Both the gate and bias voltage dependence are important. Quantum 
fluctuations due to resonant tunneling yield zero-bias anomalies as function of the bias voltage, which can be changed 
from maxima to minima by varying the gate voltage. We, furthermore, discussed the splitting of the zero-bias anomaly 
by an external magnetic field and found good agreement with recent experiments. 

We have presented a real-time approach which is based on a nonperturbative calculation of transition rates between 
different states of a local strongly correlated system coupled to fermionic or bosonic baths. We present systematic 
rules of how to set up well defined perturbation expansions for the rates in terms of the tunneling matrix elements 
between dot and leads. The formally exact rates are used to calculate occupation probabilities and the current 
from Master equations and current formulas which are intuitively obvious. The method has a wide applicability, 
ranging from the study of arbitrary dot level structures up to the investigation of macroscopic quantum coherence 
phenomena. The latter can arise from the time evolution of non-stationary initial states or by the application of 
explicitly time-dependent fields. 

The usage of real-time methods to understand low-temperature behavior of strongly correlated fermions either in 
equilibrium or non-equilibrium situations is a rather new field and has not yet been applied extensively. Compared 
with the conventional methods in imaginary timeEj they offer the possibility to set up new approximation schemes. 
In this paper we have performed a nonperturbative resummation of higher order coherent tunneling processes to 
calculate transition and current rates analytically for temperatures smaller than the intrinsic broadening T. Although 
the criterium for considering certain diagrams is yet not motivated by the usage of a "small" parameter, the diagrams 
are selected in a systematic way. We have chosen all diagrams which keep the total density matrix as close as possible 
to the diagonal state up to one electron-jhple pair excitation in the reservoirs. This reminds of techniques applied 
within variational wave function ansatzesEJ but here formulated on the basis of density matrices for nonequilibrium 
systems at finite temperatures. Furthermore, there are many possibilities to improve our approximation by considering 
more diagrams by analytical or numerical methods. Simple limiting cases as e.g. the noninteracting case are already 
exactly incorporated within our approximation. Since the strongly interacting case gives also at least qualitatively 
good results, our method may be a good candidate to cover the whole range from weak to strong interaction within 
the same approximation scheme. 

We like to thank D. Averin, J. von Delft and M. Hettler for useful discussions. The support by the Deutsche 
Forschungsgemeinschaft, through SFB 195, by the Swiss National Science Foundation (H.S.), and by the A. v. Humboldt 
award of the Academy of Finland (G.S.) is gratefully acknowledged. 
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APPENDIX A: THE RULES IN TIME SPACE 



Each term of the expansion Eq. p7| ) with operators Ai of the form Eq. ( jig ) can be calculated according to the 
following rules: 

1. Draw all topological different diagrams with directed tunneling lines connecting pairs of internal or external 

vertices containing lead electron operators. Assign a reservoir index a and a spin index a to each of these lines. 
Connect all vertices containing boson operators in all possible ways by boson lines. Assign states \ and the 
corresponding energy e x to each element of the Keldysh contour connecting two vertices. 

2. The propagation from t' to t with t' < t on the Keldysh contour implies a factor exp[— ie x (t — t')\. 

3. The state \ which is assigned to the left most part of the diagram implies a factor P° from the initial density 

matrix. Each vertex containing a dot operator B gives rise to a matrix element {x'\B\x) where x (x') i s the dot 
state entering (leaving) the vertex with respect to the Keldysh contour. 

4. Each directed tunneling line with index a running from t' to t implies (-l) v j£ (t,t') with v being the number of 

electron operators (due to external vertices) on the part of the Keldysh contour from t' to t. The line corresponds 
to a tunneling process in reservoir a. Each boson line connecting vertices at times t and t' implies P K (t,t') if 
the phase factors at these vertices have different sign. Otherwise, the boson line has the value P K (t^t 1 )^ 1 . 

5. Each diagram carries a prefactor (— i) m {— l) c , where m is the total number of internal vertices and c the number 

of crossings of tunneling lines. There may be a further minus sign due to the order of dot electron operators 
which emerges from the matrix elements (x'\B\x) discussed in rule 3. 

6. Integrate over the internal times along the Keldysh contour without changing their ordering and sum over the 

reservoir and spin indices. 

We emphasize that these diagrammatic rules hold for arbitrary dot Hamiltonians Hp — e x\x){x\i i- e -j the states 
X can be many-body eigenfunctions of Hp containing complicated correlations due to Coulomb interaction, magnetic 
fields, geometric setups, etc. . Such eigenfunctions have been calculated for special situationse3'E2l and can be used as 
an input for our diagrammatic language. In this paper, however, we will concentrate ourselves on the dot Hamiltonian 
(||) where the states x are trivially known. For this special case, the matrix elements (x'l-^lx) from rule 3 can only 
give rise to minus signs whereas they can have a more pronounced influence in more general situationscSEl 

Furthermore, we note that the same diagrammatic rules even hold for arbitrary time-dependent dot Hamiltonians 
-ffo(i) which are not diagonal in the states x- I n this case one has to assign two states x' an d x to the beginning and 
the end of each element of the Keldysh contour, respectively. The factor exp[— ie x (t — t')] from rule 2 is then replaced 
by the matrix element (x\UD(t,t')\x') where Ud denotes the time evolution operator of He, and t (t') are the times 
at the end (beginning) of the element of the Keldysh contour. 



APPENDIX B: THE RULES IN ENERGY SPACE 

We obtain the diagrammatic rules in energy space by expanding the expectation value in Eq. ( ^i"| ) and then 
performing the time integrals. We order the times of all internal (m) and external vertices (n) from left to right and 
label them by Tj with j = 1,2, ...,m + n (with r m+n = 0), irrespective on which branch they are. The Keldysh 
contour integrals are now written as ordinary integrals. This includes a minus sign for each internal vertex on the 
backward propagator. If the initial density matrix is diagonal we then encounter expressions of the type 

dn [ dT 2 ... [ dTm+n _ ie O+r le -iAS 1 (r 1 -r 2 ) e -iAS 2 (r 2 -r 3 )... e -iAB m+ „_ 1 r m+ „_ 1 



T m+ „-2 

1 1 1 



_ -m+n—1 



AE-l + i0+ AE 2 + i0+ AE m+n _x + i0+ 



(Bl) 



Here AEj is the difference of all energies going to the left minus all energies going to the right in each segment limited 
by Tj and Tj+i- This includes the energies of the propagators and - if present - the energies of the tunneling, boson 

and virtual lines. The convergence factor e 0+ri is related to an adiabatic switching on of the tunneling term Ht- The 
factor j m +™~ 1 cancels with the factor (— i) m from rule 4 above together with the prefactor (- 



Eq. (Bl). Therefore, the corresponding rules in energy representation read: 
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V . Draw all topologically different diagrams with fixed ordering of the vertices along the real axis, i.e., irrespective 
on which branch they are. The vertices are connected by tunneling and boson lines as in time space. In addition 
to the energy e x assigned to the propagators we assign an energy E to each tunneling line. For each boson line 
choose a direction (arbitrarily) and assign also an energy E. The external vertices are connected by virtual lines 
with energies Ei (i = 2, . . . , n) as described above. 

2'. For each segment derived from tj < r < Tj+\ with j = 1, 2, . . . , m + n — 1 assign a resolvent & E 1 +i0 + where AEj 
is the difference of the leftgoing minus the rightgoing energies (including the energies of the tunneling, boson 
and virtual lines). 

3'. See rule 3 in time space. 

4.' For each coupling of vertices write (—1) V ^(E), if the tunneling line of reservoir a is going backward and 
(— l) v j~(E), if it is going forward with respect to the closed time path (definition of v see rule 4 in time space). 
For each boson line write P + (E), if it is going backward and P~(E), if it is going forward with respect to the 
closed time path. 

5 ' . The prefactor is given by(— l) b (— l) c , where b is the total number of internal vertices on the backward propagator 
and c the number of crossings of tunneling lines. There may be a further minus sign due to the order of dot 
electron operators which emerges from the matrix elements discussed in rule 3. 

6'. Integrate over the energies of tunneling and boson lines and sum over the reservoir and spin indices. 



APPENDIX C: CURRENT CONSERVATION 

In this appendix we prove Eq. ( |34]) and current conservation. Let us consider any diagram E^f x (i',i) in the 
expression 

EEfW'*)- (ci) 

X V 

By changing the vertical position of the rightmost vertex we get a new diagram which has up to a minus sign the 
same value as the old diagram from which the new one was constructed. If the rightmost tunneling of the old diagram 
line has a reservoir index different from a, then the new diagram is of the form £ f ,,, so that the sum of all these 



contributions in Eq. (CI) is zero. The other diagrams are divided into two classes: in the one (the other) class, 
the rightmost tunneling line of each diagram enters (leaves) the forward propagator or leaves (enters) the backward 
propagator. The change of the vertical position of the rightmost vertex, then, increases (decreases) the value of p by 
one, so that the new diagram is of the form S^f ' }, . Furthermore, the old and the new diagram belong to different 
classes. After changing the position of the rightmost vertex of only one class and then shifting the summation index 
p in Eq. (|Cl|), we obtain exactly all diagrams of E x ^-y + x wn i cn P roves Eq. ( p4|) . 

The conservation of probability follows directly from the master equation (|2q). Summation over \ together with 
E x *V, x (t',t)=0 yields 



E 



P x (t)=0. (C2) 



To prove current or charge conservation we first recognize that 

E E x4 = - E K-,x ( C3 ) 

JV( x )=p N(x)=P+l 

where N(x) is the particle number on the dot for state x- This relation follows directly by changing the vertical 
position of the rightmost vertex. 

A fter multiplication of the master equation (|2^) with — e and N(x) and summation over x> we use Eqs. ( |33"| ) and 
( |C3|) , insert the current formula (|32| ) and find the conservation law for the total charge flowing into the dot 

52i a (t) = j t Q(t) (C4) 



1G 



where Q = —eN = —e J2 X N(x)Px 1S the charge on the dot. 

In the stationary and time-independent case Eq. (C4j) reduces to the conservation of the tunneling current 



(C5) 



whereas for the general case the r.h.s. of Eq. (C4) is minus the sum over all displacement currents flowing in the 
reservoirs. 

An importa nt re sult of this appendix is that any approximation for the rates is current conserving provided that the 
condition Eq. ( |C3| ) is satisfied. This means that we always have to consider both vertical positions of the rightmost 
vertex. 



APPENDIX D: ANALYTIC SOLUTION FOR ZERO MAGNETIC FIELD 

For zero magnetic field, i.e. e a = e for all a, we define the quantities tt(E) = n a (E), a(E) = a a (E), and 

<!>+{E) = ct>l%{a,a,E) and <t>~{E) = £ fitf (a, a', E) (Dl) 



which are independent of a. We get the integral equations 

[E-e-*(E)}<t>t(E)=± 1 ±(E)- la (E) J dE' E _^ + . Q+ ^{E') 



(D2) 



where 7 a (-B) = "f a (E) + Mj£(E) and ^(E) = ^2 a (/>a(E). Summing over a and taking the imaginary part we 
obtain the solution 



(D3) 



where we used the definitions ^(E) = £) Q 7±(£), l( E ) = + M-y+(E), 

A ± = J dE-/ ± (E)\Tr(E)\ 2 and A = / dE \tt(E)\ 2 . 

Furthermore, we obtain directly from ( |D2| ) a relation between cj) a and (f> 

l{E)^{E) = la {E)cj>±{E) ± iriEMEhtiE) - ^(E^E)} 

Using (§2|), the current rates follow from ££+ = 2iM J dElm<p+(E) and = 2i J dElm<p-(E) . With Eqs. (|d|) 
and (D5), the result is 



(D4) 



(D5) 



£«+ = -2-KiM 



^X a + J dE\n(E)\ 2 [ 1 -(E) 1 +(E)- 1 +(E) 1 -(E)} 
K-M [ dE\ir(E)\ 2 [y-(E) 1 +(E)- 7 +(E) 7 -(E)] 



(D6) 
(D7) 



where A a = / dEj a (E)\w(E)\ 2 . 

Summing the current rates over a and using J2 a A Q = A~ + MX + = 1, we get the total transition rates (note that 



A+ , A- 

Eo.o = —2iriM—— and S CT o = 2ni- 



(D8) 



A ""■ A 

and the solution of the stationary Master equation (|2^) reads 

P st = A~ and Pf = A+ with A" + M\ + = 1 . (D9) 

The stationary current follows from @ I S J = -ie[P ( fY>^ + MPfY%+] (note that = 0), which gives as the 

final result Eq. 
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FIG. 1. The SET transistor. 
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FIG. 2. A diagram showing various tunneling processes: sequential tunneling in the left and right junctions, a term preserving 
the norm, a cotunneling process, and resonant tunneling. 
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FIG. 3. The iteration of processes for the propagator II. 
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FIG. 4. Graphical representation of the current I a through lead a. Internal vertices are not indicated. 
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FIG. 5. Graphical representation of the relation between the current and the correlation functions. Here, the line connecting 
the external vertices is a real one. Internal vertices are not indicated. 
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<£> (a,a,E) = 



%2 %2 

FIG. 6. Definition of <f>. It denotes a part of a diagram with an open tunneling line entering from the right. 
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FIG. 7. The irreducible self energy is obtained by attaching the open tunneling line of cj> and (f>* to the upper and lower 
propagator. 
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FIG. 8. Graphical representation of the self-consistent equation for <f> beginning with an empty dot state. 
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FIG. 9. Graphical representation of the self-consistent equation for <j) beginning with an occupied dot state. 
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FIG. 10. The iteration of processes for the propagator it with a tunneling line running in parallel from the right to the left. 
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FIG. 11. In our approximation, the diagram for the irreducible self-energy cr <T (i5) contains one tunneling line in addition to 
the backward running line. 
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FIG. 12. Graphical representation of G&(E). 
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FIG. 13. Graphical representation ol G^{E). 




FIG. 14. The spectral density for M = 2, T = T B = 0.005r, e = -2I\ g = 0.2, o; s = 0.25r and E c = 50r at different 
voltages. For V — there are resonances at multiples of wb, which split for finite bias voltage. Inset: spectral density for 
M = 1, T = 0.00005r, T s = 0.5I\ e = -r, V = 0, p = 0.5, uj b = 0.25r and £ c = 50r. 
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FIG. 15. The differential conductance vs. bias voltage for T — Tb = 0.005r, e = -21", lob = 0.25r and Ec = 50F. The 
curves show a maximum at zero bias and satellite maxima at multiples of lob for a finite electron-boson coupling. Inset (g = 0): 
increasing voltage leads to an overall decrease of the spectral density in the range \E\ < eV . 
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FIG. 16. The differential conductance vs. bias voltage for T — Tb = 0.025r, e = 0, uj b = 0.5r and E c = 50r. The curves 
show a minimum at zero bias and steps at multiples of lob for a finite electron-boson coupling. Inset (g = 0): increasing voltage 
leads to an overall increase of the spectral density in the range \E\ < eV. 
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FIG. 17. The differential conductance as a function of e for T = 0.125I\ eV = 15V, g = 0.3, lo b = 2.5F and E c = 250V. 
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FIG. 18. The differential conductance vs. bias voltage for T = 4.3/xel/, eo-(B = 0) = —5. 2meV , F = 3AmeV, a c = 0.33, and 
-Ec = 30meV. The circles are experimental data from Ref. 
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FIG. 19. The differential conductance vs. bias voltage for T = 4.3[ieV , e CT = —5.2meV , B — 0, T = 3AmeV, a c = 0.33, and 
Ec = SOmeV. The circles are experimental data from Ref. ^J. 
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FIG. 20. The maximal linear conductance vs. temperature for e a — —5.2meV, B — 0, V = 3AmeV, a c = 0.33, and 
i?c = SOmel 7 . The circles are experimental data from Ref. 
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